###############################################################
#### Graphs and balance tables for TAPS and Canada surveys ####
###############################################################

# Load R packages
library ("haven")
library ("survey")
library ("dplyr")


# Function to plot histograms in paper
plotHist <- function (X, y.lim, main.label) {
  r <- length (X)
  breaks <- seq(from=0.5, to=(r+0.5), by=1)
  par (mar=c(1,3,3,0), las=1)
  plot (c(0.5,max(breaks)), c(0,y.lim), main=main.label
        , type="n", axes=F, ylab="", xlab="")
  axis (2)
  mtext (side=1, line=0, at=1:r, labels.o)
  for (i in 1:r){
    segments (y0=X[i], y1=X[i], x0=breaks[i], x1=breaks[i+1])
  }
  y.pos <- sort (rep(1:r,2))
  x.pos <- sort(c(1,rep(2:r,2),r+1))
  for (l in 1:(r*2)){
    segments (y0=0, y1=X[y.pos[l]], x0=breaks[x.pos[l]], x1=breaks[x.pos[l]])
  }
  text (xy.coords (1:r, X+0.5), labels=round (X,1))
  abline (h=0)
}

# Set labels to make graphs
labels.o <- c("business","politicians","global market","other")
labels.t <- c("fast/no global","slow/no global","fast/global","slow/global")



#####################
#### Canada 2015 ####
#####################

# Read data from Canada survey (Experiment 3)
setwd ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/")
Canada <- read_dta ("R&P_Canada_2015.dta")

summary (with (Canada, rowSums (cbind (jensen_a_in, jensen_b_in, jensen_c_in, jensn_d_in), na.rm=T)))

# Recode "gender"
Canada$gender <- car::recode (Canada$gender, "2=1; 1=0")
# Recode treatment category
Canada$treatment <- ifelse (!is.na(Canada$jensen_a_in), 1
                    , ifelse (!is.na(Canada$jensen_b_in), 2
                    , ifelse (!is.na(Canada$jensen_c_in), 3, 4)))

# Recode outcome of interest
Canada$creditblame <- ifelse (Canada$treatment==1, Canada$jensen_a1, 
                      ifelse (Canada$treatment==2, Canada$jensen_b1,
                      ifelse (Canada$treatment==3, Canada$jensen_c1,
                      Canada$jensen_d1)))
Canada$creditblame2 <- car::recode (Canada$creditblame,
                                    "5=NA")
Canada$creditblame2 <- as.factor (Canada$creditblame2)

# Convert 'Other' and 'Non-response' into NAs
Canada$creditblame3 <- car::recode (Canada$creditblame, "4:5=NA")




# Dummy variables used in original analyses
Canada$polresp <- ifelse (Canada$creditblame==2, 1, 0)
# 0:business people or global market forces or other or missing values
# 1:Canadian politicians
Canada$globeresp <- ifelse (Canada$creditblame==3, 1, 0)
# 0:business people or Canadian politicians or other or missing values
# 1:global market forces

# Use survey weights
svyd.Canada <- svydesign (id=~SurveyID, weights=~weight, data=Canada)
X <- svytable(formula = ~creditblame2 + treatment, design = svyd.Canada, Ntotal = 100)

# Produce plots
par (mfrow=c(2,2))
for (k in 1:4){
  plotHist (X=X[,k], y.lim=18, main.label=labels.t[k])  
}


# Prepare table with balance statistics
vars2balance.C <- c("gender","own_financ","vote2011","born","education"
                  ,"income","employment","racial","religion1","religion2"
                  ,"language","treatment")

smallCanada <- Canada[,is.element(colnames(Canada), vars2balance.C)]
smallCanada$treatment <- as.factor (smallCanada$treatment)
myTable.C <- smallCanada %>% group_by(treatment) %>%
   summarise_all(funs(m=mean(., na.rm=T), sd=sd(., na.rm=T)))

# To print table of balance statistics for Canada
CanadaTable <- as.data.frame (t(round (myTable.C[,-1],2)))
print (CanadaTable)  # Table B3



# Turn histogram plots into "confidence interval around percentages" plot
# This requires importing numbers from last table produced in R&P_Canada_TAPS_2015DoFile.do
# Population estimates (estimated percent that would hold X accountable), plus 95% confidence intervals
# for each treatment
treat1 <- rbind (c(.2686,.2045,.5269),c(.2263,.1658,.4765),c(.3156,.2496,.5767),c(2.28,2.137,2.565))
treat2 <- rbind (c(.1091,.3635,.5273),c(.0832,.3192,.4800),c(.1419,.4103,.5742),c(1.487,2.33,2.409))
treat3 <- rbind (c(.2836,.2109,.5055),c(.2409,.1721,.4558),c(.3305,.2558,.5551),c(2.287,2.134,2.54))
treat4 <- rbind (c(.1195,.3668,.5137),c(.0903,.3204,.4644),c(.1564,.4158,.5627),c(1.676,2.438,2.514))


pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1canada2015.pdf", h=10, w=7)
par (mfrow=c(4,1), mar=c(4,7,2,2), cex.axis=0.9, las=1)
plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Fast growth / No global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat1[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat1[2,1:3]*100, x1=treat1[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Slow growth / No global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat2[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat2[2,1:3]*100, x1=treat2[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Fast growth / Global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat3[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat3[2,1:3]*100, x1=treat3[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Slow growth / Global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat4[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat4[2,1:3]*100, x1=treat4[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")
dev.off ()




pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1canada2015-nonGlobalFrame.pdf", h=3, w=9)
par (mfrow=c(1,2), mar=c(4,7,2,2), cex.axis=0.8, las=1, cex.main=0.8)
plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat1[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat1[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat1[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat1[2,1:3]*100, x1=treat1[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat2[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat2[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat2[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat2[2,1:3]*100, x1=treat2[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")
dev.off ()


pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1canada2015-globalFrame.pdf", h=3, w=9)
par (mfrow=c(1,2), mar=c(4,7,2,2), cex.axis=0.8, las=1, cex.main=0.8)
plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat3[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat3[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat3[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat3[2,1:3]*100, x1=treat3[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat4[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat4[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat4[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat4[2,1:3]*100, x1=treat4[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")
dev.off()




pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1canada2015-bothFrames.pdf", h=3, w=12)
   layout (matrix(c(1,2,2,2,3,3,3), nrow=1))
   par (mar=c(3,2,3,2), cex.main=1.5, cex.axis=1.15)
   plot (c(0,1), c(0.7,3.3), type="n", axes=F, xlab="", ylab="")
   text (xy.coords (c(1,1,1), c(1,2,3))
         , labels=c("Business","Politicians","Global Markets")
         , pos=2, cex=1.3, font=2)
   plot (c(0,100), c(0.7,3.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
   axis (1)
   # axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
   text(xy.coords(treat1[1,1:3]*100,c(0.8,1.8,2.8))
        , labels=round (treat1[1,1:3]*100, 1))
   for (i in 1:3){
     text(xy.coords(treat1[1,i]*100,i-0.33)
          , labels=paste0 ("(",round (treat1[4,i], 2),")"), cex=0.8)
   }
   points (xy.coords(treat1[1,1:3]*100, c(0.95,1.95,2.95))
           , pch=19, col="gray")
   segments (x0=treat1[2,1:3]*100, x1=treat1[3,1:3]*100
             , y0=c(0.95,1.95,2.95), y1=c(0.95,1.95,2.95), lwd=2, col="gray")
   
   text(xy.coords(treat3[1,1:3]*100,c(1.30,2.30,3.30))
        , labels=round (treat3[1,1:3]*100, 1))
   for (i in 1:3){
     text(xy.coords(treat3[1,i]*100,i+0.15)
          , labels=paste0 ("(",round (treat3[4,i], 2),")"), cex=0.8)
   }
   points (xy.coords(treat3[1,1:3]*100, c(1.05,2.05,3.05))
           , pch=19, col="black")
   segments (x0=treat3[2,1:3]*100, x1=treat3[3,1:3]*100
             , y0=c(1.05,2.05,3.05), y1=c(1.05,2.05,3.05), lwd=2, col="black")
   
   
   plot (c(0,100), c(0.7,3.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
   axis (1)
   # axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
   text(xy.coords(treat2[1,1:3]*100,c(0.8,1.8,2.8))
        , labels=round (treat2[1,1:3]*100, 1))
   for (i in 1:3){
     text(xy.coords(treat2[1,i]*100,i-0.33)
          , labels=paste0 ("(",round (treat2[4,i], 2),")"), cex=0.8)
   }
   points (xy.coords(treat2[1,1:3]*100, c(0.95,1.95,2.95))
           , pch=19, col="gray")
   segments (x0=treat2[2,1:3]*100, x1=treat2[3,1:3]*100
             , y0=c(0.95,1.95,2.95), y1=c(0.95,1.95,2.95), lwd=2, col="gray")
   
   text(xy.coords(treat4[1,1:3]*100,c(1.30,2.30,3.30))
        , labels=round (treat4[1,1:3]*100, 1))
   for (i in 1:3){
     text(xy.coords(treat4[1,i]*100,i+0.15)
          , labels=paste0 ("(",round (treat4[4,i], 2),")"), cex=0.8)
   }
   points (xy.coords(treat4[1,1:3]*100, c(1.05,2.05,3.05))
           , pch=19, col="black")
   segments (x0=treat4[2,1:3]*100, x1=treat4[3,1:3]*100
             , y0=c(1.05,2.05,3.05), y1=c(1.05,2.05,3.05), lwd=2, col="black")
   legend ("bottomright", bty="n", pch=19, col=c("black","gray")
           , legend=c("'Global markets' frame", "No 'global markets' frame")
           , cex=1.3)
dev.off()


###########################
#### TAPS version 2015 ####
###########################

# Read dataset
USdata <- read_dta ("R&P_US_TAPS_2015.dta")

# Eliminate observations that are not in the design
USdata <- USdata[!is.na(USdata$xgroupab),]
summary (with (USdata, rowSums (cbind (ROSAS11S40, ROSAS21S40, ROSAS31S40, ROSAS41S40), na.rm=T)))

#### Codes are: ####
# -1: Refused
#  1: Business people
#  2: US politicians
#  3: global market forces
#  4: other
#  5: DK

# Recode relevant data
# Q1: Politicians are responsible
USdata$polresp <- ifelse ((USdata$ROSAS11S40==2 & !is.na(USdata$ROSAS11S40)) |
                          (USdata$ROSAS21S40==2 & !is.na(USdata$ROSAS21S40)) |
                          (USdata$ROSAS31S40==2 & !is.na(USdata$ROSAS31S40)) |
                          (USdata$ROSAS41S40==2 & !is.na(USdata$ROSAS41S40)), 1, 
                          ifelse ((USdata$ROSAS11S40==-1 & !is.na(USdata$ROSAS11S40)) |
                                  (USdata$ROSAS21S40==-1 & !is.na(USdata$ROSAS21S40)) |
                                  (USdata$ROSAS31S40==-1 & !is.na(USdata$ROSAS31S40)) |
                                  (USdata$ROSAS41S40==-1 & !is.na(USdata$ROSAS41S40)), NA, 0))


# Q2: Global forces are responsible
USdata$globeresp <- ifelse ((USdata$ROSAS11S40==3 & !is.na(USdata$ROSAS11S40)) |
                            (USdata$ROSAS21S40==3 & !is.na(USdata$ROSAS21S40)) |
                            (USdata$ROSAS31S40==3 & !is.na(USdata$ROSAS31S40)) |
                            (USdata$ROSAS41S40==3 & !is.na(USdata$ROSAS41S40)), 1, 
                          ifelse ((USdata$ROSAS11S40==-1 & !is.na(USdata$ROSAS11S40)) |
                                    (USdata$ROSAS21S40==-1 & !is.na(USdata$ROSAS21S40)) |
                                    (USdata$ROSAS31S40==-1 & !is.na(USdata$ROSAS31S40)) |
                                    (USdata$ROSAS41S40==-1 & !is.na(USdata$ROSAS41S40)), NA, 0))

# Q2: Politicians are responsible
USdata$polresp2 <- ifelse ((USdata$ROSAS13S40==2 & !is.na(USdata$ROSAS13S40)) |
                            (USdata$ROSAS23S40==2 & !is.na(USdata$ROSAS23S40)) |
                            (USdata$ROSAS33S40==2 & !is.na(USdata$ROSAS33S40)) |
                            (USdata$ROSAS43S40==2 & !is.na(USdata$ROSAS43S40)), 1, 
                          ifelse ((USdata$ROSAS13S40==-1 & !is.na(USdata$ROSAS13S40)) |
                                    (USdata$ROSAS23S40==-1 & !is.na(USdata$ROSAS23S40)) |
                                    (USdata$ROSAS33S40==-1 & !is.na(USdata$ROSAS33S40)) |
                                    (USdata$ROSAS43S40==-1 & !is.na(USdata$ROSAS43S40)), NA, 0))

# Q2: Global forces are responsible
USdata$globeresp2 <- ifelse ((USdata$ROSAS13S40==3 & !is.na(USdata$ROSAS13S40)) |
                             (USdata$ROSAS23S40==3 & !is.na(USdata$ROSAS23S40)) |
                             (USdata$ROSAS33S40==3 & !is.na(USdata$ROSAS33S40)) |
                             (USdata$ROSAS43S40==3 & !is.na(USdata$ROSAS43S40)), 1, 
                           ifelse ((USdata$ROSAS13S40==-1 & !is.na(USdata$ROSAS13S40)) |
                                     (USdata$ROSAS23S40==-1 & !is.na(USdata$ROSAS23S40)) |
                                     (USdata$ROSAS33S40==-1 & !is.na(USdata$ROSAS33S40)) |
                                     (USdata$ROSAS43S40==-1 & !is.na(USdata$ROSAS43S40)), NA, 0))

# Q2: Credit/blame
USdata$creditblame <- ifelse (USdata$xgroupab==1, USdata$ROSAS13S40, 
                              ifelse (USdata$xgroupab==2, USdata$ROSAS23S40,
                                      ifelse (USdata$xgroupab==3, USdata$ROSAS33S40,
                                              USdata$ROSAS43S40)))
USdata$creditblame <- car::recode (USdata$creditblame, "-1=NA; 5=NA")

# Generate control variables
USdata$age <- USdata$agesp
USdata$female <- ifelse (USdata$gendersp==1 & !is.na(USdata$gendersp), 0
                         , ifelse (USdata$gendersp==2 & !is.na(USdata$gendersp), 1, NA))
USdata$Obamaapproval <- car::recode (USdata$apprpressp, "-1=NA; 5=NA") #1:strongly approve; 4:strongly disapprove
USdata$democrat <- ifelse (USdata$PARTYID1SP==1 & !is.na(USdata$PARTYID1SP), 1, 0)
USdata$republican <- ifelse (USdata$PARTYID1SP==2 & !is.na(USdata$PARTYID1SP), 1, 0)
USdata$independent <- ifelse (USdata$PARTYID1SP==3 & !is.na(USdata$PARTYID1SP), 1, 0)
USdata$treatment <- as.factor (USdata$xgroupab)

#### Survey commands ####
svyd.US <- svydesign (id=~wustlid, weights=~mar2015wt1, data=USdata)
Y <- svytable(formula = ~creditblame + xgroupab, design = svyd.US, Ntotal = 100)

# Plot graph of results for TAPS 2015
par (mfrow=c(2,2))
for (k in 1:4){
   plotHist (X=Y[,k], y.lim=18, main.label=labels.t[k])  
}


#### Balance statistics ####
vars2balance <- c("age","female","Obamaapproval","democrat","republican"
                  ,"independent","incomehhldsp","maritalstatsp","religattendsp","treatment")

USdatasmall <- USdata[,is.element(colnames(USdata), vars2balance)]

myTable <- USdatasmall %>% group_by(treatment) %>%
   summarise_all(funs(m=mean(., na.rm=T), sd=sd(., na.rm=T)))

# To print table
USTable <- as.data.frame (t(round (myTable[,-1],2)))
print (USTable)  # Table B2




# Turn histogram plots into "confidence interval around percentages" plot
# This requires importing numbers from last table produced in R&P_US_TAPS_2015DoFile.do
# Population estimates (estimated percent that would hold X accountable), plus 95% confidence intervals
# for each treatment
treat1 <- rbind (c(.4182,.1697,.4121),c(.3296,.0949,.3202),c(.5124,.2848,.5107),c(4.71,4.791,4.911))
treat2 <- rbind (c(.0994,.5778,.3228),c(.0626,.4867,.2447),c(.1542,.6639,.4123),c(2.291,4.563,4.304))
treat3 <- rbind (c(.4485,.1664,.3851),c(.3628,.1153,.3040),c(.5375,.2341,.4731),c(4.497,3.013,4.348))
treat4 <- rbind (c(.2778,.3697,.3524),c(.2383,.3238,.3089),c(.3213,.4181,.3986),c(4.434,5.234,4.633))


pdf ("experiment1usa2015.pdf", h=10, w=7)
par (mfrow=c(4,1), mar=c(4,7,2,2), cex.axis=0.9, las=1)
plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Fast growth / No global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat1[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat1[2,1:3]*100, x1=treat1[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Slow growth / No global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat2[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat2[2,1:3]*100, x1=treat2[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Fast growth / Global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat3[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat3[2,1:3]*100, x1=treat3[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,70), c(0.5,3), type="n", axes=F, xlab="", ylab="", main="Slow growth / Global frame")
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
points (xy.coords(treat4[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat4[2,1:3]*100, x1=treat4[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")
dev.off ()




pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1usa2015-nonGlobalFrame.pdf", h=3, w=9)
par (mfrow=c(1,2), mar=c(4,7,2,2), cex.axis=0.8, las=1, cex.main=0.8)
plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat1[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat1[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat1[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat1[2,1:3]*100, x1=treat1[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat2[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat2[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat2[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat2[2,1:3]*100, x1=treat2[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")
dev.off ()


pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1usa2015-globalFrame.pdf", h=3, w=9)
par (mfrow=c(1,2), mar=c(4,7,2,2), cex.axis=0.8, las=1, cex.main=0.8)
plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat3[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat3[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat3[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat3[2,1:3]*100, x1=treat3[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")

plot (c(0,100), c(0.5,3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
axis (1)
axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat4[1,1:3]*100,c(0.8,1.8,2.8)), labels=round (treat4[1,1:3]*100, 1), cex=0.8)
points (xy.coords(treat4[1,1:3]*100,1:3), pch=19, col="gray")
segments (x0=treat4[2,1:3]*100, x1=treat4[3,1:3]*100, y0=1:3, y1=1:3, lwd=2, col="gray")
dev.off()


pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1usa2015-bothFrames.pdf", h=3, w=12)
layout (matrix(c(1,2,2,2,3,3,3), nrow=1))
par (mar=c(3,2,3,2), cex.main=1.5, cex.axis=1.15)
plot (c(0,1), c(0.7,3.3), type="n", axes=F, xlab="", ylab="")
text (xy.coords (c(1,1,1), c(1,2,3))
      , labels=c("Business","Politicians","Global Markets")
      , pos=2, cex=1.3, font=2)

plot (c(0,100), c(0.7,3.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
axis (1)
# axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat1[1,1:3]*100,c(0.85,1.85,2.85))
     , labels=round (treat1[1,1:3]*100, 1))
for (i in 1:3){
   text(xy.coords(treat1[1,i]*100,i-0.33)
        , labels=paste0 ("(",round (treat1[4,i], 2),")"), cex=0.8)
}
points (xy.coords(treat1[1,1:3]*100, c(0.95,1.95,2.95))
        , pch=19, col="gray")
segments (x0=treat1[2,1:3]*100, x1=treat1[3,1:3]*100
          , y0=c(0.95,1.95,2.95), y1=c(0.95,1.95,2.95), lwd=2, col="gray")

text(xy.coords(treat3[1,1:3]*100,c(1.33,2.33,3.33))
     , labels=round (treat3[1,1:3]*100, 1))
for (i in 1:3){
   text(xy.coords(treat3[1,i]*100,i+0.15)
        , labels=paste0 ("(",round (treat3[4,i], 2),")"), cex=0.8)
}
points (xy.coords(treat3[1,1:3]*100, c(1.05,2.05,3.05))
        , pch=19, col="black")
segments (x0=treat3[2,1:3]*100, x1=treat3[3,1:3]*100
          , y0=c(1.05,2.05,3.05), y1=c(1.05,2.05,3.05), lwd=2, col="black")


plot (c(0,100), c(0.7,3.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Who should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
axis (1)
# axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
text(xy.coords(treat2[1,1:3]*100,c(0.85,1.85,2.85))
     , labels=round (treat2[1,1:3]*100, 1))
for (i in 1:3){
   text(xy.coords(treat2[1,i]*100,i-0.33)
        , labels=paste0 ("(",round (treat2[4,i], 2),")"), cex=0.8)
}
points (xy.coords(treat2[1,1:3]*100, c(0.95,1.95,2.95))
        , pch=19, col="gray")
segments (x0=treat2[2,1:3]*100, x1=treat2[3,1:3]*100
          , y0=c(0.95,1.95,2.95), y1=c(0.95,1.95,2.95), lwd=2, col="gray")

text(xy.coords(treat4[1,1:3]*100,c(1.33,2.33,3.33))
     , labels=round (treat4[1,1:3]*100, 1))
for (i in 1:3){
   text(xy.coords(treat4[1,i]*100,i+0.15)
        , labels=paste0 ("(",round (treat4[4,i], 2),")"), cex=0.8)
}
points (xy.coords(treat4[1,1:3]*100, c(1.05,2.05,3.05))
        , pch=19, col="black")
segments (x0=treat4[2,1:3]*100, x1=treat4[3,1:3]*100
          , y0=c(1.05,2.05,3.05), y1=c(1.05,2.05,3.05), lwd=2, col="black")
legend ("bottomright", bty="n", pch=19, col=c("black","gray")
        , legend=c("'Global markets' frame", "No 'global markets' frame")
        , cex=1.3)
dev.off()

###########################
#### TAPS version 2014 ####
###########################
USdata2 <- read_dta ("R&P_US_TAPS_2014.dta")

# Recode treatment and outcomes
USdata2$treatment <- as.factor (USdata2$xrosas)

USdata2$credit1 <- ifelse ((USdata2$ROSAS1S29==1 & !is.na(USdata2$ROSAS1S29)) |
                              (USdata2$ROSAS1S29==2 & !is.na(USdata2$ROSAS1S29)), 1,
                               ifelse ((USdata2$ROSAS1S29==3 & !is.na(USdata2$ROSAS1S29)) |
                                          (USdata2$ROSAS1S29==4 & !is.na(USdata2$ROSAS1S29)), 0, NA))

USdata2$blame2 <- ifelse ((USdata2$ROSAS11S29==1 & !is.na(USdata2$ROSAS11S29)) |
                              (USdata2$ROSAS11S29==2 & !is.na(USdata2$ROSAS11S29)), 1,
                           ifelse ((USdata2$ROSAS11S29==3 & !is.na(USdata2$ROSAS11S29)) |
                                      (USdata2$ROSAS11S29==4 & !is.na(USdata2$ROSAS11S29)), 0, NA))

USdata2$credit3 <- ifelse ((USdata2$ROSAS21S29==1 & !is.na(USdata2$ROSAS21S29)) |
                              (USdata2$ROSAS21S29==2 & !is.na(USdata2$ROSAS21S29)), 1,
                           ifelse ((USdata2$ROSAS21S29==3 & !is.na(USdata2$ROSAS21S29)) |
                                      (USdata2$ROSAS21S29==4 & !is.na(USdata2$ROSAS21S29)), 0, NA))

USdata2$blame4 <- ifelse ((USdata2$ROSAS31S29==1 & !is.na(USdata2$ROSAS31S29)) |
                              (USdata2$ROSAS31S29==2 & !is.na(USdata2$ROSAS31S29)), 1,
                           ifelse ((USdata2$ROSAS31S29==3 & !is.na(USdata2$ROSAS31S29)) |
                                      (USdata2$ROSAS31S29==4 & !is.na(USdata2$ROSAS31S29)), 0, NA))

USdata2$creditblame5 <- ifelse ((USdata2$ROSAS41S29==1 & !is.na(USdata2$ROSAS41S29)) |
                             (USdata2$ROSAS41S29==2 & !is.na(USdata2$ROSAS41S29)), 1,
                          ifelse ((USdata2$ROSAS41S29==3 & !is.na(USdata2$ROSAS41S29)) |
                                     (USdata2$ROSAS41S29==4 & !is.na(USdata2$ROSAS41S29)), 0, NA))

USdata2$creditblame <- ifelse (is.na(USdata2$credit1) &
                               is.na(USdata2$blame2) &
                               is.na(USdata2$credit3) &
                               is.na(USdata2$blame4) &
                               is.na(USdata2$creditblame5), NA, 
                               ifelse ((USdata2$credit1==1 & !is.na(USdata2$credit1)) |
                                          (USdata2$blame2==1 & !is.na(USdata2$blame2)) |
                                          (USdata2$credit3==1 & !is.na(USdata2$credit3)) |
                                          (USdata2$blame4==1 & !is.na(USdata2$blame4)) |
                                          (USdata2$creditblame5==1 & !is.na(USdata2$creditblame5)), 1, 0)) 

# Retrospective credit, ordinal
USdata2$retcreditblame <- ifelse (!is.na(USdata2$credit1), USdata2$credit1,
                                  ifelse (!is.na(USdata2$blame2), USdata2$blame2,
                                          ifelse (!is.na(USdata2$credit3), USdata2$credit3,
                                                  ifelse (!is.na(USdata2$blame4), USdata2$blame4,
                                                          ifelse (!is.na(USdata2$creditblame5), USdata2$creditblame5, NA)))))

USdata2$retcreditblameOrd <- ifelse (!is.na(USdata2$ROSAS1S29), USdata2$ROSAS1S29,
                                  ifelse (!is.na(USdata2$ROSAS11S29), USdata2$ROSAS11S29,
                                          ifelse (!is.na(USdata2$ROSAS21S29), USdata2$ROSAS21S29,
                                                  ifelse (!is.na(USdata2$ROSAS31S29), USdata2$ROSAS31S29,
                                                          ifelse (!is.na(USdata2$ROSAS41S29), USdata2$ROSAS41S29, NA)))))

USdata2$retcreditblameOrd <- car::recode (USdata2$retcreditblameOrd, "-1=NA")

# Dummy for great deal of credit or blame on retro 
USdata2$credit1g <- ifelse ((USdata2$ROSAS1S29==1 & !is.na(USdata2$ROSAS1S29)), 1,
                            ifelse ((USdata2$ROSAS1S29==2 & !is.na(USdata2$ROSAS1S29)) |
                                       (USdata2$ROSAS1S29==3 & !is.na(USdata2$ROSAS1S29))|
                                       (USdata2$ROSAS1S29==4 & !is.na(USdata2$ROSAS1S29)), 0, NA))

USdata2$blame2g <- ifelse ((USdata2$ROSAS11S29==1 & !is.na(USdata2$ROSAS11S29)), 1,
                            ifelse ((USdata2$ROSAS11S29==2 & !is.na(USdata2$ROSAS11S29)) |
                                       (USdata2$ROSAS11S29==3 & !is.na(USdata2$ROSAS11S29))|
                                       (USdata2$ROSAS11S29==4 & !is.na(USdata2$ROSAS11S29)), 0, NA))

USdata2$credit3g <- ifelse ((USdata2$ROSAS21S29==1 & !is.na(USdata2$ROSAS21S29)), 1,
                           ifelse ((USdata2$ROSAS21S29==2 & !is.na(USdata2$ROSAS21S29)) |
                                      (USdata2$ROSAS21S29==3 & !is.na(USdata2$ROSAS21S29))|
                                      (USdata2$ROSAS21S29==4 & !is.na(USdata2$ROSAS21S29)), 0, NA))

USdata2$blame4g <- ifelse ((USdata2$ROSAS31S29==1 & !is.na(USdata2$ROSAS31S29)), 1,
                           ifelse ((USdata2$ROSAS31S29==2 & !is.na(USdata2$ROSAS31S29)) |
                                      (USdata2$ROSAS31S29==3 & !is.na(USdata2$ROSAS31S29))|
                                      (USdata2$ROSAS31S29==4 & !is.na(USdata2$ROSAS31S29)), 0, NA))

USdata2$creditblame5g <- ifelse ((USdata2$ROSAS41S29==1 & !is.na(USdata2$ROSAS41S29)), 1,
                            ifelse ((USdata2$ROSAS41S29==2 & !is.na(USdata2$ROSAS41S29)) |
                                       (USdata2$ROSAS41S29==3 & !is.na(USdata2$ROSAS41S29))|
                                       (USdata2$ROSAS41S29==4 & !is.na(USdata2$ROSAS41S29)), 0, NA))

# generate variable for credit or blame for great deal of credit or blame */
USdata2$retcreditblameg <- ifelse (!is.na(USdata2$credit1g), USdata2$credit1g,
                                  ifelse (!is.na(USdata2$blame2g), USdata2$blame2g,
                                          ifelse (!is.na(USdata2$credit3g), USdata2$credit3g,
                                                  ifelse (!is.na(USdata2$blame4g), USdata2$blame4g,
                                                          ifelse (!is.na(USdata2$creditblame5g), USdata2$creditblame5g, NA)))))

# Dummy for any credit or blame on prospective credit or blame */
USdata2$procredit1 <- ifelse ((USdata2$ROSAS3S29==1 & !is.na(USdata2$ROSAS3S29)) |
                                 (USdata2$ROSAS3S29==2 & !is.na(USdata2$ROSAS3S29)) , 1,
                            ifelse ((USdata2$ROSAS3S29==3 & !is.na(USdata2$ROSAS3S29)) |
                                       (USdata2$ROSAS3S29==4 & !is.na(USdata2$ROSAS3S29)), 0, NA))

USdata2$problame2 <- ifelse ((USdata2$ROSAS13S29==1 & !is.na(USdata2$ROSAS13S29)) |
                                 (USdata2$ROSAS13S29==2 & !is.na(USdata2$ROSAS13S29)) , 1,
                              ifelse ((USdata2$ROSAS13S29==3 & !is.na(USdata2$ROSAS13S29)) |
                                         (USdata2$ROSAS13S29==4 & !is.na(USdata2$ROSAS13S29)), 0, NA))

USdata2$procredit3 <- ifelse ((USdata2$ROSAS23S29==1 & !is.na(USdata2$ROSAS23S29)) |
                                (USdata2$ROSAS23S29==2 & !is.na(USdata2$ROSAS23S29)) , 1,
                             ifelse ((USdata2$ROSAS23S29==3 & !is.na(USdata2$ROSAS23S29)) |
                                        (USdata2$ROSAS23S29==4 & !is.na(USdata2$ROSAS23S29)), 0, NA))


USdata2$problame4 <- ifelse ((USdata2$ROSAS33S29==1 & !is.na(USdata2$ROSAS33S29)) |
                                (USdata2$ROSAS33S29==2 & !is.na(USdata2$ROSAS33S29)) , 1,
                             ifelse ((USdata2$ROSAS33S29==3 & !is.na(USdata2$ROSAS33S29)) |
                                        (USdata2$ROSAS33S29==4 & !is.na(USdata2$ROSAS33S29)), 0, NA))

# Dummy for great deal of credit or blame on prospective credit or blame */
USdata2$procredit1g <- ifelse ((USdata2$ROSAS3S29==1 & !is.na(USdata2$ROSAS3S29)), 1,
                            ifelse ((USdata2$ROSAS3S29==2 & !is.na(USdata2$ROSAS3S29)) |
                                       (USdata2$ROSAS3S29==3 & !is.na(USdata2$ROSAS3S29))|
                                       (USdata2$ROSAS3S29==4 & !is.na(USdata2$ROSAS3S29)), 0, NA))

USdata2$problame2g <- ifelse ((USdata2$ROSAS13S29==1 & !is.na(USdata2$ROSAS13S29)), 1,
                               ifelse ((USdata2$ROSAS13S29==2 & !is.na(USdata2$ROSAS13S29)) |
                                          (USdata2$ROSAS13S29==3 & !is.na(USdata2$ROSAS13S29))|
                                          (USdata2$ROSAS13S29==4 & !is.na(USdata2$ROSAS13S29)), 0, NA))

USdata2$procredit3g <- ifelse ((USdata2$ROSAS23S29==1 & !is.na(USdata2$ROSAS23S29)), 1,
                               ifelse ((USdata2$ROSAS23S29==2 & !is.na(USdata2$ROSAS23S29)) |
                                          (USdata2$ROSAS23S29==3 & !is.na(USdata2$ROSAS23S29))|
                                          (USdata2$ROSAS23S29==4 & !is.na(USdata2$ROSAS23S29)), 0, NA))

USdata2$problame4g <- ifelse ((USdata2$ROSAS33S29==1 & !is.na(USdata2$ROSAS33S29)), 1,
                              ifelse ((USdata2$ROSAS33S29==2 & !is.na(USdata2$ROSAS33S29)) |
                                         (USdata2$ROSAS33S29==3 & !is.na(USdata2$ROSAS33S29))|
                                         (USdata2$ROSAS33S29==4 & !is.na(USdata2$ROSAS33S29)), 0, NA))

USdata2$procreditblame5g <- ifelse ((USdata2$ROSAS43S29==1 & !is.na(USdata2$ROSAS43S29)), 1,
                              ifelse ((USdata2$ROSAS43S29==2 & !is.na(USdata2$ROSAS43S29)) |
                                         (USdata2$ROSAS43S29==3 & !is.na(USdata2$ROSAS43S29))|
                                         (USdata2$ROSAS43S29==4 & !is.na(USdata2$ROSAS43S29)), 0, NA))

USdata2$procreditblameg <- ifelse (!is.na(USdata2$procredit1g), USdata2$procredit1g,
                                   ifelse (!is.na(USdata2$problame2g), USdata2$problame2g,
                                           ifelse (!is.na(USdata2$procredit3g), USdata2$procredit3g,
                                                   ifelse (!is.na(USdata2$problame4g), USdata2$problame4g,
                                                           ifelse (!is.na(USdata2$procreditblame5g), USdata2$procreditblame5g, NA)))))

# Prospective vote intentions 
USdata2$votepro <- ifelse ((USdata2$ROSAS4S29 != -1 & !is.na(USdata2$ROSAS4S29)), USdata2$ROSAS4S29,
                           ifelse ((USdata2$ROSAS14S29 != -1 & !is.na(USdata2$ROSAS14S29)), USdata2$ROSAS14S29,
                           ifelse ((USdata2$ROSAS24S29 != -1 & !is.na(USdata2$ROSAS24S29)), USdata2$ROSAS24S29,
                           ifelse ((USdata2$ROSAS34S29 != -1 & !is.na(USdata2$ROSAS34S29)), USdata2$ROSAS34S29,
                           ifelse ((USdata2$ROSAS43S29 != -1 & !is.na(USdata2$ROSAS43S29)), USdata2$ROSAS43S29, NA)))))

USdata2$votepro2 <- ifelse (!is.na(USdata2$votepro), 0, NA)
USdata2$votepro2 <- ifelse ((USdata2$votepro==1 & !is.na(USdata2$votepro)) |
                               (USdata2$votepro==2 & !is.na(USdata2$votepro)), 1, USdata2$votepro2)

# Generating Control Variables 
USdata2$dem <- ifelse (!is.na(USdata2$PARTYID1S27), 0, NA)
USdata2$dem <- ifelse ((USdata2$PARTYID1S27==1 & !is.na(USdata2$PARTYID1S27)), 1, USdata2$dem)

USdata2$rep <- ifelse (!is.na(USdata2$PARTYID1S27), 0, NA)
USdata2$rep <- ifelse ((USdata2$PARTYID1S27==2 & !is.na(USdata2$PARTYID1S27)), 1, USdata2$rep)

USdata2$ind <- ifelse (!is.na(USdata2$PARTYID1S27), 0, NA)
USdata2$ind <- ifelse ((USdata2$PARTYID1S27==3 & !is.na(USdata2$PARTYID1S27)), 1, USdata2$ind)

USdata2$globalization <- ifelse ((!is.na(USdata2$treatment) & as.character(USdata2$treatment)==3) |
                                   (!is.na(USdata2$treatment) & as.character(USdata2$treatment)==4), 1, 
                                 ifelse (is.na(USdata2$treatment), NA, 0)) 


#### Survey commands ####
USdata2 <- USdata2[!is.na(USdata2$treatment),]

svyd.US2 <- svydesign (id=~wustlid, weights=~apr2014wt1, data=USdata2)
Y <- svytable(formula = ~retcreditblameOrd + treatment, design = svyd.US2, Ntotal = 100)

# Slightly altered function to plot results from US 2014
newPlotHist <- function (X, y.lim, main.label, axislabel) {
   r <- length (X)
   breaks <- seq(from=0.5, to=(r+0.5), by=1)
   par (mar=c(3,2,3,0), las=1)
   plot (c(0.5,max(breaks)), c(0,y.lim), main=main.label
         , type="n", axes=F, ylab="", xlab="")
   axis (2)
   mtext (side=1, line=0, at=1:r, axislabel, cex=0.7)
   mtext (side=1, line=1, at=c(1,4,5), c("more","less","DK"), cex=0.7)
   mtext (side=1, line=2, at=2.5, c("<---- responsibility ---->"), cex=0.7)
   for (i in 1:r){
      segments (y0=X[i], y1=X[i], x0=breaks[i], x1=breaks[i+1])
   }
   y.pos <- sort (rep(1:r,2))
   x.pos <- sort(c(1,rep(2:r,2),r+1))
   for (l in 1:(r*2)){
      segments (y0=0, y1=X[y.pos[l]], x0=breaks[x.pos[l]], x1=breaks[x.pos[l]])
   }
   text (xy.coords (1:r, X+0.5), labels=round (X,1))
   abline (h=0)
}

# Slight change of labels to correspond to 2014
labels.2014 <- c(labels.t, "control")
labels.x <- 1:5

# Make plot, Figure C1
par (mfrow=c(1,5))
for (k in 1:5){
   newPlotHist (X=Y[,k], y.lim=10, main.label=labels.2014[k], axislabel=labels.x)  
}


#### Balance statistics US 2014 ####
USdata2$age <- USdata2$agesp
USdata2$female <- ifelse (USdata2$gendersp==1 & !is.na(USdata2$gendersp), 0
                         , ifelse (USdata2$gendersp==2 & !is.na(USdata2$gendersp), 1, NA))
USdata2$Obamaapproval <- car::recode (USdata2$apprpressp, "-1=NA; 5=NA") #1:strongly approve; 4:strongly disapprove
USdata2$democrat <- ifelse (USdata2$PARTYID1SP==1 & !is.na(USdata2$PARTYID1SP), 1, 0)
USdata2$republican <- ifelse (USdata2$PARTYID1SP==2 & !is.na(USdata2$PARTYID1SP), 1, 0)
USdata2$independent <- ifelse (USdata2$PARTYID1SP==3 & !is.na(USdata2$PARTYID1SP), 1, 0)

vars2balance <- c("age","female","Obamaapproval","democrat","republican"
                  ,"independent","incomehhldsp","maritalstatsp"
                  ,"religattendsp","treatment")

USdata2small <- USdata2[,is.element(colnames(USdata2), vars2balance)]

myTable2 <- USdata2small %>% group_by(treatment) %>%
   summarise_all(funs(m=mean(., na.rm=T), sd=sd(., na.rm=T)))

# To print table
USTable2 <- as.data.frame (t(round (myTable2[,-1],2)))
print (USTable2)  # Table B1



# Turn histogram plots into "confidence interval around percentages" plot
# This requires importing numbers from bottom rows of last tables produced in R&P_US_TAPS_2014DoFile.do
# Population estimates (estimated percent that would hold X accountable), plus 95% confidence intervals
# for each treatment
retro1 <- c(2.62,1.291,5.246,0.9382)
retro2 <- c(47.54,36.97,58.33,5.526)
retro3 <- c(2.544,1.182,5.392,0.9864)
retro4 <- c(50.28,39.32,61.21,5.67)
retro5 <- c(44.66,34.91,54.84,5.146)

prosp1 <- c(2.857,1.337,5.998,1.095)
prosp2 <- c(39.89,30.07,50.6,5.302)
prosp3 <- c(3.927,1.923,7.855,1.413)
prosp4 <- c(48.36,37.39,59.49,5.725)
prosp5 <- c(39.68,28.8,51.68,5.93)

incumb1 <- c(37.75,28.77,47.66,4.867)
incumb2 <- c(27.95,20.26,37.21,4.348)
incumb3 <- c(31.19,21.89,42.3,5.258)
incumb4 <- c(30.14,21.68,40.2,4.761)
incumb5 <- c(36.08,27.63,45.48,4.593)


pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1usa2014-retrospective.pdf", h=3, w=12)
# layout (matrix(c(1,2,2,2,3,3,3), nrow=1))
par (mfrow=c(1,2), mar=c(3,3,3,2), cex.main=1.5, cex.axis=1.15)
# plot (c(0,1), c(0.7,2.3), type="n", axes=F, xlab="", ylab="")
# text (xy.coords (c(1,1), c(1,2))
#       , labels=c("Retrospective","Prospective")
#       , pos=2, cex=1.3, font=2)
plot (c(0,100), c(0.7,1.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Politicians should be", bold (" credited "), "for", bold (" fast "), "growth?", sep=" ")))
axis (1)
# axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
# retro domestic fast
text(xy.coords(retro1[1],c(0.85))
     , labels=round (retro1[1], 1))
text(xy.coords(retro1[1],c(0.78))
     , labels=paste0("(",round (retro1[4], 3),")"), cex=0.8)
points (xy.coords(retro1[1], c(0.9))
        , pch=19, col="gray")
segments (x0=retro1[2], x1=retro1[3]
          , y0=c(0.9), y1=c(0.9), lwd=2, col="gray")

# retro global fast
text(xy.coords(retro3[1],c(1.25))
     , labels=round (retro3[1], 1))
text(xy.coords(retro3[1],c(1.17))
     , labels=paste0("(",round (retro3[4], 3),")"), cex=0.8)
points (xy.coords(retro3[1], c(1.1))
        , pch=19, col="black")
segments (x0=retro3[2], x1=retro3[3]
          , y0=c(1.1), y1=c(1.1), lwd=2, col="black")

# retro control
# text(xy.coords(retro5[1],c(0.95))
#      , labels=round (retro5[1], 1))
# points (xy.coords(retro5[1], c(1))
#         , pch=4, col="black")
# segments (x0=retro5[2], x1=retro5[3]
#           , y0=c(1), y1=c(1), lty=2, col="black")

# # prosp domestic fast
# text(xy.coords(prosp1[1],c(1.8))
#      , labels=round (prosp1[1], 1))
# points (xy.coords(prosp1[1], c(1.9))
#         , pch=19, col="gray")
# segments (x0=prosp1[2], x1=prosp1[3]
#           , y0=c(1.9), y1=c(1.9), lwd=2, col="gray")
# 
# # prosp global fast
# text(xy.coords(prosp3[1],c(2.2))
#      , labels=round (prosp3[1], 1))
# points (xy.coords(prosp3[1], c(2.1))
#         , pch=19, col="black")
# segments (x0=prosp3[2], x1=prosp3[3]
#           , y0=c(2.1), y1=c(2.1), lwd=2, col="black")
# 
# # prosp control
# text(xy.coords(prosp5[3]+2,c(2))
#      , labels=round (prosp5[1], 1), pos=2)
# points (xy.coords(prosp5[1], c(2))
#         , pch=4, col="black")
# segments (x0=prosp5[2], x1=prosp5[3]
#           , y0=c(2), y1=c(2), lty=2, col="black")

plot (c(0,100), c(0.7,1.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Politicians should be", bold (" blamed "), "for", bold (" slow "), "growth?", sep=" ")))
axis (1)
# axis (2, at=c(1,2,3), tick=FALSE, labels=c("Business","Politicians","Global Markets"), cex=0.8)
# retro domestic slow
text(xy.coords(retro2[1],c(0.85))
     , labels=round (retro2[1], 1))
text(xy.coords(retro2[1],c(0.78))
     , labels=paste0("(",round (retro2[4], 3),")"), cex=0.8)
points (xy.coords(retro2[1], c(0.9))
        , pch=19, col="gray")
segments (x0=retro2[2], x1=retro2[3]
          , y0=c(0.9), y1=c(0.9), lwd=2, col="gray")

# retro domestic fast
text(xy.coords(retro4[1],c(1.25))
     , labels=round (retro4[1], 1))
text(xy.coords(retro4[1],c(1.17))
     , labels=paste0("(",round (retro4[4], 3),")"), cex=0.8)
points (xy.coords(retro4[1], c(1.1))
        , pch=19, col="black")
segments (x0=retro4[2], x1=retro4[3]
          , y0=c(1.1), y1=c(1.1), lwd=2, col="black")

# retro control
# text(xy.coords(retro5[1]+2,c(0.95))
#      , labels=round (retro5[1], 1))
# points (xy.coords(retro5[1], c(1))
#         , pch=4, col="black")
# segments (x0=retro5[2], x1=retro5[3]
#           , y0=c(1), y1=c(1), lty=2, col="black")

# # prosp domestic slow
# text(xy.coords(prosp2[1],c(1.8))
#      , labels=round (prosp2[1], 1))
# points (xy.coords(prosp2[1], c(1.9))
#         , pch=19, col="gray")
# segments (x0=prosp2[2], x1=prosp2[3]
#           , y0=c(1.9), y1=c(1.9), lwd=2, col="gray")
# 
# # prosp domestic fast
# text(xy.coords(prosp4[1],c(2.2))
#      , labels=round (prosp4[1], 1))
# points (xy.coords(prosp4[1], c(2.1))
#         , pch=19, col="black")
# segments (x0=prosp4[2], x1=prosp4[3]
#           , y0=c(2.1), y1=c(2.1), lwd=2, col="black")
# 
# # prosp control
# text(xy.coords(prosp5[3]+2,c(2))
#      , labels=round (prosp5[1], 1), pos=2)
# points (xy.coords(prosp5[1], c(2))
#         , pch=4, col="black")
# segments (x0=prosp5[2], x1=prosp5[3]
#           , y0=c(2), y1=c(2), lty=2, col="black")

dev.off()



pdf ("~/Dropbox/Openness_Survival/ReplicationPackageRandP/experiment1usa2014-incumbent.pdf", h=3, w=12)
par (mfrow=c(1,2), mar=c(3,2,3,2), cex.main=1.5, cex.axis=1.15)
plot (c(0,100), c(0.7,1.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Plans to vote for incumbent under", bold (" fast "), "growth", sep=" ")))
axis (1)
# incumb domestic fast
text(xy.coords(incumb1[1],c(0.85))
     , labels=round (incumb1[1], 1))
text(xy.coords(incumb1[1],c(0.78))
     , labels=paste0 ("(", round (incumb1[4], 1), ")"), cex=0.8)
points (xy.coords(incumb1[1], c(0.9))
        , pch=19, col="gray")
segments (x0=incumb1[2], x1=incumb1[3]
          , y0=c(0.9), y1=c(0.9), lwd=2, col="gray")
# incumb global fast
text(xy.coords(incumb3[1],c(1.25))
     , labels=round (incumb3[1], 1))
text(xy.coords(incumb3[1],c(1.17))
     , labels=paste0 ("(", round (incumb3[4], 1), ")"), cex=0.8)
points (xy.coords(incumb3[1], c(1.1))
        , pch=19, col="black")
segments (x0=incumb3[2], x1=incumb3[3]
          , y0=c(1.1), y1=c(1.1), lwd=2, col="black")
# incumb control
# text(xy.coords(incumb5[1],c(0.95))
#      , labels=round (incumb5[1], 1))
# points (xy.coords(incumb5[1], c(1))
#         , pch=4, col="black")
# segments (x0=incumb5[2], x1=incumb5[3]
#           , y0=c(1), y1=c(1), lty=2, col="black")

plot (c(0,100), c(0.7,1.3), type="n", axes=F, xlab="", ylab="", main=expression(paste("Plans to vote for incumbent under", bold (" slow "), "growth?", sep=" ")))
axis (1)
# incumb domestic slow
text(xy.coords(incumb2[1],c(0.85))
     , labels=round (incumb2[1], 1))
text(xy.coords(incumb2[1],c(0.78))
     , labels=paste0 ("(", round (incumb2[4], 1), ")"), cex=0.8)
points (xy.coords(incumb2[1], c(0.9))
        , pch=19, col="gray")
segments (x0=incumb2[2], x1=incumb2[3]
          , y0=c(0.9), y1=c(0.9), lwd=2, col="gray")
# incumb domestic fast
text(xy.coords(incumb4[1],c(1.25))
     , labels=round (incumb4[1], 1))
text(xy.coords(incumb4[1],c(1.17))
     , labels=paste0 ("(", round (incumb4[4], 1), ")"), cex=0.8)
points (xy.coords(incumb4[1], c(1.1))
        , pch=19, col="black")
segments (x0=incumb4[2], x1=incumb4[3]
          , y0=c(1.1), y1=c(1.1), lwd=2, col="black")
# incumb control
# text(xy.coords(incumb5[1]+2,c(0.95))
#      , labels=round (incumb5[1], 1))
# points (xy.coords(incumb5[1], c(1))
#         , pch=4, col="black")
# segments (x0=incumb5[2], x1=incumb5[3]
#           , y0=c(1), y1=c(1), lty=2, col="black")
# legend ("bottomright", bty="n", pch=c(19,19,4), col=c("black","gray",
#                                               "black")
#         , legend=c("'Global markets' frame", "'Domestic' frame", "Control")
#         , cex=1.3)
legend ("bottomright", bty="n", pch=c(19,19), col=c("black","gray")
        , legend=c("'Global markets' frame", "'Domestic' frame")
        , cex=1.3)
dev.off()

